﻿#pragma once

#include <learn/alg/alg_geom.h>
#include <learn/geom/geom_bspline_curve.h>
#include <learn/geom/geom_bspline_surface.h>

namespace xi
{

namespace alg
{

using namespace xi::geom;

template <std::size_t N, bool Rat = false> struct Alg_Surface_Common
{
    using Vec = VecT<N, double>;        // P
    using VecR = VecT<N + Rat, double>; // P'

    /**
     * @brief 构造线性插值角点的样条曲面
     *
     * @param[in] c0 S(0,0)
     * @param[in] c1 S(1,0)
     * @param[in] c2 S(0,1)
     * @param[in] c3 S(1,1)
     * @return std::shared_ptr<Geom_BSplineSurface<N, Rat>>
     */
    static std::shared_ptr<Geom_BSplineSurface<N, Rat>> bilinear(const Vec &c0, const Vec &c1, const Vec &c2,
                                                                 const Vec &c3)
    {
        std::vector<double> uKnots = {0, 0, 1, 1};
        std::vector<double> vKnots = {0, 0, 1, 1};

        if constexpr (Rat)
        {
            std::vector<VecR> cpts = {homogenous_to_cartesian(c0), homogenous_to_cartesian(c3),
                                      homogenous_to_cartesian(c1), homogenous_to_cartesian(c2)};
            return std::make_shared<Geom_BSplineSurface<N, Rat>>(std::move(cpts), std::move(uKnots), std::move(vKnots),
                                                                 1, 1);
        }
        else
        {
            std::vector<VecR> cpts = {c0, c3, c1, c2};
            return std::make_shared<Geom_BSplineSurface<N, Rat>>(std::move(cpts), std::move(uKnots), std::move(vKnots),
                                                                 1, 1);
        }
    }

    /**
     * @brief 样条曲线沿着指定方向平移得到圆柱
     *
     * @param[in] C
     * @param[in] d
     * @return std::shared_ptr<Geom_BSplineSurface<N, Rat>>
     */
    static std::shared_ptr<Geom_BSplineSurface<N, Rat>> cylinder(std::shared_ptr<Geom_BSplineCurve<N, Rat>> C,
                                                                 const Vec &d)
    {
        auto &ccpts = C->poles();
        std::vector<double> uKnots = {0.0, 0.0, 1.0, 1.0};
        std::vector<double> vKnots = C->knots();
        std::vector<VecR> cpts(C->pole_num() * 2);
        for (std::size_t i = 0; i < C->pole_num(); i++)
        {
            auto w = C->weight(i);
            cpts[i] = cpts[i + C->pole_num()] = ccpts[i];
            for (std::size_t j = 0; j < N; j++)
                cpts[i][j] += d[j] * w;
        }
        return std::make_shared<Geom_BSplineSurface<N, Rat>>(std::move(cpts), std::move(uKnots), std::move(vKnots), 1,
                                                             C->degree(), false);
    }

    /**
     * @brief 计算两条样条曲线的直纹面
     *
     * @param[in] C
     * @param[in] D
     * @return std::shared_ptr<Geom_BSplineSurface<N, Rat>>
     */
    static std::shared_ptr<Geom_BSplineSurface<N, Rat>> ruled(std::shared_ptr<Geom_BSplineCurve<N, Rat>> C,
                                                              std::shared_ptr<Geom_BSplineCurve<N, Rat>> D)
    {
        Alg_Knot<N>::compat<Rat>(C, D);

        auto &ccpts = C->poles();
        auto &dcpts = D->poles();
        std::vector<double> uKnots = {0.0, 0.0, 1.0, 1.0};
        std::vector<double> vKnots = C->knots();
        std::vector<VecR> cpts(C->pole_num() * 2);
        for (std::size_t i = 0; i < C->pole_num(); i++)
        {
            cpts[i] = ccpts[i];
            cpts[i + C->pole_num()] = dcpts[i];
        }
        return std::make_shared<Geom_BSplineSurface<N, Rat>>(std::move(cpts), std::move(uKnots), std::move(vKnots), 1,
                                                             C->degree());
    }

    /**
     * @brief 一条样条曲线绕一个轴旋转 theta 得到的样条曲面
     *
     * @param[in] C 初始曲线
     * @param[in] S 原点
     * @param[in] T 轴向
     * @param[in] theta 转角
     * @return std::shared_ptr<Geom_BSplineSurface<N, true>>
     */
    static std::shared_ptr<Geom_BSplineSurface<N, true>> revolution(std::shared_ptr<Geom_BSplineCurve<N, Rat>> C,
                                                                    const Vec &S, const Vec &T, double theta)
    {
        std::vector<double> uKnots(12);
        auto pi_2 = math::Numeric::pi / 2;
        std::size_t narcs = 1;
        theta = degrees_to_radians(theta);
        if (theta > pi_2 && theta <= pi_2 * 2)
        {
            narcs = 2;
            uKnots[3] = uKnots[4] = 0.5;
        }
        else if (theta <= pi_2 * 3)
        {
            narcs = 3;
            uKnots[3] = uKnots[4] = 1.0 / 3;
            uKnots[5] = uKnots[6] = 2.0 / 3;
        }
        else
        {
            narcs = 4;
            uKnots[3] = uKnots[4] = 0.25;
            uKnots[5] = uKnots[6] = 0.5;
            uKnots[7] = uKnots[8] = 0.75;
        }

        auto dtheta = theta / narcs;
        auto j = 3 + 2 * (narcs - 1);
        for (std::size_t i = 0; i < 3; j++, i++)
        {
            uKnots[i] = 0;
            uKnots[j] = 1;
        }
        uKnots.resize(j);

        auto wm = std::cos(dtheta / 2);
        auto angle = 0.0;

        std::vector<double> cosines(narcs + 1);
        std::vector<double> sines(narcs + 1);
        for (std::size_t i = 1; i < narcs + 1; i++)
        {
            angle += dtheta;
            cosines[i] = std::cos(angle);
            sines[i] = std::sin(angle);
        }

        auto m = C->pole_num() - 1;
        std::vector<VecT<N + 1, double>> cpts((uKnots.size() - 3) * (m + 1));
        for (std::size_t j = 0; j < m + 1; j++)
        {
            auto P0 = C->pole(j);
            auto O = project_to_line(P0, S, T);
            auto X = P0 - O;
            auto Y = cross(T, X);
            auto T0 = Y;

            auto index = 0;
            auto w = C->weight(j);
            cpts[j] = cartesian_to_homogenous(P0, w);
            for (std::size_t i = 1; i < narcs + 1; i++)
            {
                auto P2 = O + X * cosines[i] + Y * sines[i];
                cpts[(index + 2) * C->pole_num() + j] = cartesian_to_homogenous(P2, w);

                auto T2 = Y * cosines[i] - X * sines[i];
                auto option = intersect_line_line(P0, T0, P2, T2);
                if (!option.has_value())
                {
                    // 此时说明 P0, P2 重合
                    cpts[(index + 1) * C->pole_num() + j] = cartesian_to_homogenous(P0, wm * w);
                }
                else
                {
                    auto [alpha, beta] = option.value();
                    cpts[(index + 1) * C->pole_num() + j] = cartesian_to_homogenous(P0 + T0 * alpha, wm * w);
                }

                if (i < narcs)
                {
                    P0 = P2;
                    T0 = T2;
                }
                index += 2;
            }
        }

        std::vector<double> vKnots = C->knots();
        return std::make_shared<Geom_BSplineSurface<N, true>>(std::move(cpts), std::move(uKnots), std::move(vKnots), 2,
                                                              C->degree());
    }

    static std::shared_ptr<Geom_BSplineSurface<N>> spherical()
    {
        return nullptr;
    }
};

} // namespace alg

} // namespace xi